Linear response formula for finite frequency thermal conductance of open systems 
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An exact linear response expression is obtained for the heat current in a classical Hamiltonian 
system coupled to heat baths with time-dependent temperatures. The expression is equally valid 
at zero and finite frequencies. We present numerical results on the frequency dependence of the 
response function for three different one-dimensional models of coupled oscillators connected to 
Langevin baths with oscillating temperatures. For momentum conserving systems, a low frequency 
peak is seen that, is higher than the zero frequency response for large systems. For momentum 
non-conserving systems, there is no low frequency peak. The momentum non-conserving system 
is expected to satisfy Fourier's law, however, at the single bond level, we do not see any clear 
agreement with the predictions of the diffusion equation even at low frequencies. We also derive an 
exact analytical expression for the response of a chain of harmonic oscillators to a (not necessarily 
small) temperature difference; the agreement with the linear response simulation results for the 
same system is excellent. 

PACS numbers: 



I. INTRODUCTION 

In many low dimensional systems, heat transport unex- 
pectedly violates Fourier's law of heat conduction [ll-Q- 
This can be because of integrability or proximity to in- 
tegr ability, which is more common in low dimensions, as 
recognized starting from the Fermi-Pasta-Ulam (FPU) 
model Q. Alternatively, even ergodic low-dimensional 
systems can show anomalous heat conduction, with the 
conductivity diverging with system size, if they conserve 
momentum. Apart from the theoretical interest, under- 
standing heat transport in such systems is of relevance 
to heat conduction in carbon nanotubes 

Most of the recent activity @, Q in this field has dealt 
with the zero-frequency conductivity. But time depen- 
dent temperature sources have been discussed in exper- 
imental situations in the context of measuring the fre- 
quency dependent thermal conductivity @, 0] and spe- 
cific heat [8( of glassy systems. Theoretically, there have 
been a few studies on the frequency dependent thermal 
current response using a microscopic approach based on 
Luttinger's derivation of the Green-Kubo formula and a 
hypothesis about the equality of certain transport coef- 
ficients Q, and from a phenomenological approach [To| . 
A recent paper studied thermal ratchet effects in an in- 
homogeneous anharmonic chain coup led to baths with 
time-dependent temperatures [TJ G3 - 

In this paper, we adopt a different approach: we find 
the linear heat conductance of a system placed in contact 
with two heat reservoirs with time-dependent tempera- 
tures. Physically the notion of bath temperatures oscil- 
lating in time make sense if we assume that the frequency 
of oscillation is much smaller compared to time scales for 
local thermal equilibration in the reservoirs. An exact 
expression (in the linear response regime) for the heat 
current due to a small oscillating temperature difference 



between the reservoirs is obtained. 

Our earlier result [l3j obtained the zero frequency con- 
ductance of a finite system rather than the conductivity 
in the infinite system limit. Thus the thermodynamic 
limit was not taken first (in fact, not at all), in contrast to 
the standard Green-Kubo formula [T3|, which cannot be 
applied when the infinite system conductivity diverges. 
Our expression for the zero frequency conductance in- 
volved the heat current auto-correlation function for an 
open system. The extension to finite frequencies in this 
paper follows the same approach, with the response now 
depending on the position inside the system where the 
current is measured. 

We also show results of numerical simulations for the 
frequency dependent response function by measuring the 
appropriate correlation function. For one-dimensional 
momentum conserving anharmonic crystals, we find a 
resonant response at a frequency uj ~ 1/JV for a chain 
of N particles due to sound waves propagating from one 
end of the system to the other. As N increases, the reso- 
nance gets broader and its height decreases slightly. How- 
ever, its height relative to the zero-frequency response in- 
creases, and for large N this resonance is stronger than 
the zero frequency response. 

We find that the low frequency peak disappears for sys- 
tems where momentum is not conserved. Fourier's law 
is known to be valid for such systems, so that the heat 
current should satisfy the diffusion equation. If one com- 
pares the numerical results for the frequency-dependent 
heat current with the prediction from the diffusion equa- 
tion at the single bond level, there seem to be substantial 
discrepancies. 

Numerical simulations for the frequency dependent re- 
sponse function of a one-dimensional harmonic crystal, 
and an exact analytical expression for the full response 
(for finite AT) of the same, are also presented. As far as 
we are aware of this is the first example of a case where an 
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analytical expression for the response function has been 
obtained. For a harmonic system the full response is also 
linear and hence we expect the linear response result to 
agree with the exact response function. Indeed we find 
excellent agreement between the numerical simulations of 
the expression of the linear response and the numerically 
evaluated exact response expression. 

All three systems mentioned above also show a high- 
frequency peak in the response function, whose location 
is independent of N. One can loosely ascribe this to the 
fact that the dynamics in the interior of the system are 
underdamped (actually, undamped), so that particles ap- 
proaching each other recoil, and the heat current auto- 
correlation function shows rapid oscillations in the tem- 
poral domain. Such high-frequency oscillations are not 
seen in hard particle models, such as the Random Col- 
lision Model [l5|]. This is discussed further when we de- 
rive the analytical expression for the harmonic oscilla- 
tor. However, a quantitative understanding of the high- 
frequency peak is lacking. 



II. OSCILLATOR CHAINS WITH LANGEVIN 
BATHS 

We follow the derivation of Ref. [l3| to obtain the fi- 
nite frequency heat conductance of an oscillator chain 
with Langevin baths at the ends; more detail is provided 
in Ref. [13]. Consider the motion of N particles on a one 
dimensional lattice, described by the following Hamilto- 
nian: 
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where x = {xi} and v = {vi} with I — 1,2, ...N are 
the displacements of the particles about their equilibrium 
positions and their velocities, and {mi} are their masses. 
We assume fixed boundary conditions, xq — xn+i — 
0. The particles 1 and N are connected to white noise 
Langevin heat baths at temperatures Tj, and Tr. Thus 
the equations of motion are 
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_d_ 

dxi 



[U(xi-i - xi) + U(xi - xi+i) + V(xi)} 



+ Si t i[r]L(t) - 7lwi] + Si yN [r] R (t) - ^rvn] (2) 



for I = 1,2, ...N. Here T]L t n(t) are uncorrelated zero 
mean Gaussian noise terms satisfying the fluctuation dis- 
sipation relations 

{VL, R {t)VLM(t'))v = 27L,Rk B T LiR S(t - t') , (3) 

where denotes an average over the noise. 

The derivation of the linear response theory starts with 
the Fokker-Planck equation for the full phase space dis- 
tribution function P(x;v;i). If T L = T R = T, the 
steady state solution to the equation is the equilibrium 



Boltzmann distribution. We now assume that the tem- 
peratures at the two ends are oscillating in time with 
Tl,r = T ± AT(i)/2. We will obtain a perturbative 
solution about the equilibrium solution. The steps are 
very similar to the standard derivation of the fluctuation 
dissipation theorem. The Fokker Planck equation corre- 
sponding to Eq. (J2j) is 



dxi 



OiP + O n P 
(4) 

where /; = —dH/dxi is the force acting on the Vth parti- 
cle. The operators Oi.n come from the Langevin damp- 
ing and noise on the terminal particles: 
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N 
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(5) 



With Tl,r = T±AT(t)/2, we can group terms according 
to their power of AT to obtain 



dP 
~dt 
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L AT P 



where 
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For AT = 0, the steady state solution of the Fokker 
Planck equation is the equilibrium Boltzmann distribu- 
tion Pq = exp[— f3H]/Z, where Z is the canonical parti- 
tion function and (3 = l/(fesT). For AT =/= 0, we start 
with the equilibrium distribution at time t = to and then 
let the system evolve under the full Fokker Planck oper- 
ator. Writing P(x, v, t) = Pq + p(x, v, t) and retaining 
terms to O(AT), 



dp 

at 



= Lp + L AT P { 



(8) 



Setting to - 
tion 

p(x;v;t) 



-oo we get the formal solution to this equa- 



,{t-t')l A/3(t') J /p (v)P (x,v)dt' (9) 



where J/ P (v) is defined by 

OP 
~dt 
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(A/3)J /P P 



(10) 



P=Pn 



from which 

1r 



J fp 



1L 



2mN 



[mNV N -k B T]-- [miv 1 ~k B T]. (11) 
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The expectation value of any function (A A) = (A) 
{A)o of any observable ^4(x; v) then takes the form: 



(AA(t)) AT = - 



1 



k B T 2 



{A( T )J fp (0)}AT(t - r)dr 



(12) 



3 



where we have defined the equilibrium average 
(A(t)J fp (0)) = fdxf dvAe Lt J fp P and we have used 
the time translational invariance of the equilibrium cor- 
relation function. In particular, we are interested in the 
energy current between two adjacent particles . The in- 
stantaneous current from the Vth to the I + l'th site 
is given by: = \{vi + where = 

—dU(xi — xi+i) / dxi+i is the force on the I + l'th parti- 
cle due to the l'th particle. We get for the average heat 
current flowing between any bond on the chain by: 



1 



AT 



k B T 2 Jo 



{jl+X,l(T)J fp (0))AT(t-T)dT. 

(13) 



For a oscillating temperature given by AT(t) 
AT(w)e luJt this gives: 



AT(uj)e iujt 
1 



{ji+iArUfpm^dr , (14) 



where Gi(u) is the magnitude of the response — to be 
computed numerically in Section IIIII — and cj>i is the 
phase. The correlation function (ji+i,i(T)Jf p (Q)) on the 
right hand side of this equation is for a system in equi- 
librium at temperature T. 

A few comments are appropriate here. First, as shown 
in Ref . [l3| , for u> — > it is possible to manipulate the in- 
tegrand on the right and make it proportional to the aut- 
correlation function of the heat current integrated over 
the entire chain, J^i 3l+i,l( T )i yielding a result resembling 
the standard Green-Kubo formula (but without the ther- 
modynamic limit). This manipulation is not possible for 
u/0. Thus the current response depends on /, the po- 
sition inside the chain where the response is measured, 
as one would expect. Moreover, the correlation function 
involves J/ p , which is different from the heat current. 

Second, although we have assumed that ATl = — ATr 
to resemble the zero-frequency calculations of Ref. [l3| 
where such an assumption is appropriate, at u ^ there 
is no reason why one cannot treat AT^ and ATr as in- 
dependent variables. It is straightforward to extend the 
derivation above and obtain the response to ATr and 
ATl , with Jf p in Eq. (|T4"]) replaced by the first and second 
part of Eq. (fTT|) respectively. For large N, one expects 
that the response to a oscillatory temperature perturba- 
tion at one end of the chain should only depend on the 
distance from that end and be the same as for a semi- 
infinite chain. 

Finally, expressions similar to Eq. (fT3"|) can be obtained 
for any quantity that depends on the phase space vari- 
ables of the system, not just j;+ij(r). It does not apply to 
the heat current flowing into the system from the reser- 
voirs, since they involve the Langevin noise T)l,r, and 
these have to be obtained indirectly. Thus Eq. (fTB")) is 



valid for I = 1, and one also has 
1 d 



(&i(t)/dt> 



AT 



k B T 2 dt Jq 



ei(T)J /p (0)>AT(t-r)dr. 

Replacing the dj dt with a —d/dr acting on AT and inte- 
grating by parts, adding this to Eq. (H"3l) . and using the 
fact that j2i(i) + de\{t)/dt = ji,L{t) (where jx t L is the 
heat current flowing in from the left reservoir), we have 
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AT = 



k B T 2 Jo 
1 



k B T 2 



{ji,L(r)Jf P (0))AT(t-T)dr 
AT(t)( ei (0)J /p (0)). (16) 



Fourier transforming, for AT(t) = AT(uj)e lu}t , the heat 
current flowing from the left reservoir is 



■7i,lM 



1 



k B T 2 



(jix(r)Jf P (0)) 



z-^ T dT+^k B . 



(17) 

This response function has a non-zero uj — > oo limit 
from the second term on the right hand side. This is 
reasonable: if AT oscillates at a very high frequency, 
the effect on (x, v) should be negligible, but the current 
flowing from the left reservoir should oscillate because 
(VL(t)vi(t)) v = ^hksT^it^jmx is proportional to the in- 
stantaneous temperature of the reservoir. The instanta- 
neous response of Eq. (fTTf is a peculiarity of white noise 
stochastic baths, and is not seen for Nose-Hoover baths 
- where even the heat current at the boundary is in 
terms of the extended phase space variables — or a fluid 
system with Maxwell boundary conditions where conti- 
nuity requires that the heat current at the boundary and 
just inside the system should be the same. Therefore, 
hereafter we work with j'21 and Jn,n-i when we want 
the current at the boundaries. 

Although the derivation given above is for a one- 
dimensional chain, it is straightforward to see that it is 
valid for any system that is connected to only two reser- 
voirs, regardless of its dimensionality. 



III. NUMERICAL RESULTS 

Numerical simulations to obtain the correlation func- 
tion on the right hand side of Eq. (TT4|) were performed 
on three different systems, which differ in the potential 
of each particle. From these correlation functions we ob- 
tained Gi(lu) using Eq. (JUJ). The velocity- Ver let algo- 
rithm was used, with a time step St = 0.005. We veri- 
fied that doubling St does not change our results. For 
the largest systems, the initial equilibration time was 
t eq = 64 x 10 6 , after which the dynamical equations were 
evolved for a time t = 5 x 10 s . All the particle masses 
were set to 1, 71, = 7# = 1, and the reservoirs were at 
temperature T = 2.0. Figured] shows Gi(ui) as a function 
of w, as defined by Eq. (jl"4")l . for FPU chains of different 
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FIG. 1: (Color online) Plot of magnitude of the response func- 
tion, Gi(uj), for FPU chains of different lengths. The inset 
shows the correlation function Ci(t), which has the same in- 
formation in the time domain. 



FIG. 2: (Color online) Plot of the magnitude of the response 
function, G2(w), for FPU-chains of different lengths. The 
inset shows Gi(uu) for various N = 64 and various I. 



lengths. The potential used was U(x) = x 2 /2 + x 4 /4 
with V(x) = 0. An N- independent high frequency peak 
and a low frequency peak at uj ~ 1/N are seen. Higher 
harmonics of the low frequency peak can be barely dis- 
cerned. As the system size is increased, the low fre- 
quency peak broadens and decreases slightly in height, 
but the zero frequency response drops much faster. Thus 
by N — 128, the uj ~ 1/A resonance is clearly stronger 
than the zero frequency response. Note that Eq. (fT4|) 
gives the conductance, not the conductivity; the ui = 
conductance decreases as ~ l/A 1- ". It is expected that 
a = 1/3 Q for large N but this would require much 
larger system sizes to verify. The inset to Figure [T] shows 
CiW = (j2i(t)Jf p (0)), i.e. the same information in the 
time domain. N- independent short time oscillations that 
decay to (approximately) zero are seen. An 'echo' of the 
oscillation is seen at a time r^r that is approximately 
N/v, where v is possibly related to the velocity of effec- 
tive phonons [l6| . At high frequencies, G\{uj) is approx- 
imately independent of N as one would expect, with a 
high frequency peak. As w — > oo, Gi(w) ~ 1/uS 2 . 

Figure [5] shows G2(w), the magnitude of the response 
function at a distance I = 2 from the left boundary. The 
low frequency peak (and its harmonics) are still present, 
but much more irregular in shape. However, from a de- 
vice perspective, it is the currents flowing into the bound- 
aries that are important. The high frequency behavior is 
independent of N, and as seen in the inset, the peak in 
Gi(ui) shifts to smaller uj as I is increased. It is not clear 
if G2(w ~ oo) ~ 1/uj 6 as is seen for the harmonic chain 
(discussed later in this paper). 

Figure [3] shows G\ (ui) for chains of different lengths 
with an onsite potential V(x) = x 4 /4. The interparticle 
potential is harmonic, U(x) = x 2 /2. The dynamics are 
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FIG. 3: (Color online) Plot of the magnitude of the response 
function, Gi(cj), for <^> 4 -chains of different lengths. The inset 
shows Ci(t). 



not momentum conserving, and the zero frequency con- 
ductance should be inversely proportional to N. This is 
not seen in the data for two reasons: direct measurement 
of the zero frequency conductance by applying a small 
temperature difference between the reservoirs shows that 
one needs N > 256 to see the ~ 1/N dependence, and 
the curves for the two larger systems (more noticeably 
N = 128) have not reached their w — >• limit in the 
figure. The low frequency resonance is gone, replaced 
by a broad A-independent plateau. This is presumably 
because at finite temperature, the effective phonons are 
optical instead of acoustic. The A-independent high fre- 
quency peak is also present. The response in the interior 
of the chain, shown in Figure U is similar, except that 
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FIG. 4: (Color onine) Plot of Gi(u) for </> 4 -chains of different 
lengths. A fit to ~ 1/lj 6 in the asymptotic high frequency 
regime is shown. The inset has Gi(ui) for various I and N = 
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FIG. 5: (Color online) Plot of Gi(w) for <fi 4 chains of different 
lengths (LR) and Gi iB (uj) from the diffusion equation (DE). 



the low frequency plateau extends down to uj — (or to 
very small uj). As for the FPU chains, we fit G\(uj ~ oo 
to ~ 1/uj 2 and — less successfully — G2(w — > oo) to 
~ 1/ui 6 . From the inset to Figure 01 there is no signifi- 
cant Z-dependence to the location of the high frequency 
peak in Gi(u), unlike what we saw for FPU chains. 

Beyond the ~ l/N dependence of the zero frequency 
conductance, one expects that heat transport in systems 
that are not momentum conserving should be diffusive, 
and the temperature field will satisfy dTi/dt = D(Ti + \ — 
2Ti + where D = k/C is the diffusion constant. 

With an ~ e lut time dependence, the resultant differ- 
ence equation can be solved with T^(uj) and Tr(u>) spec- 
ified, and thence the heat current ji+i,z = k(Ti — Xj+i) 



FIG. 6: (Color online) Plot of Gi(w) for harmonic chains 
of different lengths, from the analytical expression derived in 
Section llVl Because of the complicated structure in the figure, 
N = 128 is not included. The linear response simulation 
results for N = 64 are also shown (LR). The inset shows 
Ci(t). 



can be calculated. Some features of the solution are 
Gf iS (oj = 0) cx l/N, Gf iS {uj) is independent of N for 
N -> oo, Gf iS {uj -> 0) - w 1 / 2 exp[-(w/2L>) 1 / 2 Z] and 
Gf lS (oj — > oo) ~ 1/uj 1 . In Figure [5] we plot the responses 
Gi lS together with the linear response results G\ for the 
<jr model. For each system size we fix the diffusion con- 
stant D so that the ui = results for the two responses 
match. One expects that the low-frequency agreement 
between the two sets should become better with increas- 
ing system size. However this is not clear from our data. 
At high frequencies, the expectation Gi*(ui) ~ l/u> is 
definitely not borne out. Since the diffusion equation is 
not expected to be valid at microscopic time or length 
scales, and the fact that ~ l/N scaling of the zero fre- 
quency heat conductance is only seen for N > 256 sug- 
gests that 'microscopic' length scales are quite large here, 
the lack of agreement at the single bond level and high 
frequencies is perhaps not surprising. A clear under- 
standing of this requires further work. 

Finally, we show the results for a harmonic chain, with 
V(x) = and U(x) = x 2 /2. In this case we show in the 
next section [sec. (|IV|) ] that the response can be 

obtained exactly and expressed in terms of a single inte- 
gral over frequencies. Here we give numerical results for 
Gi(ui) obtained using this exact formula [Eq. (|21l) ] and 
also compare it with the linear response result [Eq. (IT4|) ] . 
We show G\{ui) in Figure |51 with results from numeri- 
cal simulations of the linear response formula also shown 
for N — 64. We see excellent agreement between the 
analytical and linear response result. One can see that 
Gi(u) = 0) is almost A^-independent as expected, and the 
low frequency resonance and its harmonics are more pro- 
nounced than for the FPU chain, which is not surprising 
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FIG. 7: (Color online) Plot of (?2(w) and Ga{oj) for harmonic 
chains of various lengths. The fit to the asymptotic form 
G;(u — > oo) ~ 1/u; 41 ' 2 is shown. The inset shows Gi(u>) for 
various I and N = 64. 



since there is no dispersion or damping in the interior of 
the chain. The high frequency peak seems to be present 
but is difficult to cleanly separate from the low frequency 
structure. As was the case for the FPU and </3 4 chains, 
Gi(uj — ► oo) is N- independent and ~ 1/uj 2 . In Section llVl 
the asymptotic form Gi(uj oo) ~ l/w 4 '~ 2 is obtained. 
Figure [7] shows G 2 (ui) for various system sizes, with all 
features as expected. 



IV. RESPONSE OF A HARMONIC CHAIN 

Although the integrability of the harmonic oscillator 
chain makes its behavior non-generic, and its applica- 
bility to physical systems limited, the advantage of this 
model is that its response can be completely obtained 
analytically (with some integrals evaluated numerically) 
and compared to the simulation results. We now proceed 
with the analysis. 

In this case both V(x) and U(x) are quadratic and 
the Hamiltonian can be written in the form H = 
XMX/2 + X&X/2 where M and $ are respectively 
the mass matrix and the force-constant matrix for the 
system. We will obtain the solution of the equations 
of motion in the time-dependent steady state by us- 
ing Fourier transforms in the time domain. The ap- 
proach is similar to that used in the derivation of the 
Landauer-type formula for steady state heat current in 
harmonic systems, where the current is expressed in 
terms of phonon Green's functions [17]. Let us intro- 
duce the transforms: xi(Q) = (l/2w) dtxi{t)e lUt and 

VlAW = (V 271 ") JZc dtrj LtR (t)e int . Then the Fourier 
transform solution of Eqns. © gives: 



where G + (Sl) = [-Mil 2 + $ - E+(fi)]- 1 is the phonon 
Green's function [TtJ and £ + , the self-energy correction 
due to baths, is a N x N matrix whose only non-zero 
elements are = iil^L, E^-jv = *^7-R- The noise cor- 
relations corresponding to the oscillating temperatures 
T L =T + AT/ 2 cos wi, T R = T - AT/ 2 coswi are given 
by: 

(fj L (ni)fj L (Q 2 )) = ^{ T (5(^1 + fla) 

+ (at/4)[ 5(0! + a 2 + u) + <y(«i + n 2 - u) } } , 

(rte(fii)*te(0 2 )) = ^{ Tgfa + n 2 ) 

- (AT/4)[ <5(fii + n 2 +uj) + 5(Sli + n 2 - w) ] } (19) 

and f/L , fj R are uncorrelated. The heat current on any 
bond is given by the noise average (ji+i.i) — ((l/2)(fc(x/ — 
xi+i)(vi + vi+i)), where k is the force constant of the 
bonds, and thus involves evaluating 



/oo />oo 
dfti / dfl 2 (-<n 2 ) e 
-oo J —oo 



(20) 



and this is readily evaluated using the noise properties in 
Eq. (p~9|) . After some simplifications we finally obtain: 



1 f°° 

GAuj) = \ — dnn 

x [7£{ei|i(n-w)-^ M (n-a;)} 
x{g+(-0) + g+ +11 (-f7)} 

- -/ R {gl N (fi-uj) -g+ +1N (n-u)} 



.(21) 



For nearest neighbor interactions, the force matrix $ is a 
tri-diagonal matrix. Using the properties of inverse of a 
tri-diagonal matrix we can explicitly evaluate the Green's 
function elements that are required. For simplicity con- 
sider the case k = 1 and 7l = J R = 7- Let us define A^ m 
as the determinant of the sub- matrix of [— Mf2 2 + <£> — E + ] 
that starts from the / th row and column and ends in the 
TO th row and column. We also define D^ m as the deter- 
minant of the sub-matrix of [—Mil 2 + $] starting from 
the Z th row and column and ending in the m th row and 
column. In terms of these one has: 



A, 



+ l,N 



Ai. 



/v 



i-i 



Ar. 



TV 



with 



xi(si) = g+(n)fj L (si) + g+ N (n)fj R (si) , 



(18) 



Ai,;_i = Di t i-i - iVL^D 2 ,i-\ 
Ai+i,jv = A+i,jv _ ^7^+i,jv-i 
A x ,jv = D hN - i07(£>i,JV-i + -D 2 ,at) - n 2 <y 2 D 2 , N -i ■ 

For an ordered harmonic chain with all masses equal to 
1 it is easy to show that -D;. m = sin(m — I + 2)q/ sinq 



7 



where ft 2 — 2(1 — cosg). Using this it is easy to numer- 
ically evaluate the response function Gi(ui) in Eq. (|2"Tj) 
for given values of I, N . We show some numerical results 
in Figs. (|6I7[) where we have also compared with results 
from simulations for the linear response. As expected 
the exact response and the linear response give almost 
identical results. However we have not been able to ana- 
lytically show the equivalence of the exact response and 
the linear response expressions. 

For large f2, we have q ~ n + i In £1 2 , hence 



2\l 



!/(-&) 



(22) 



This can also be seen from the equations of motion: when 
D, » 0, the dynamical equations become —miVl 2 xi = 
kxi-i. The boundary condition is — miQ 2 Xi = 771, (O), 
in which the right hand sign is effectively unity when 
calculating the Green's function. Combining these equa- 
tions, we obtain Eq. (|22|) . But a ~ l/(—Q 2 ) 1 dependence 
at large frequencies implies that the 2Z'th derivative of 
Q\ i At) has a (^-function at the origin, i.e. QtAf) ~ t 21 ^ 1 
for t > (This can be verified directly in the time do- 
main: Xi(t) oc t 21 ^ 1 satisfies the equations of motion for 
t > 0.) But then in the time domain, Eq. (f2~Tj) is equiv- 
alent to Gi(t) oc Q^' 1 (t)dtG[~i(t) for t > 0, where we have 
assumed that I is in the left half of the chain. Therefore 
Gi{t > 0) oc t 41 - 3 . Since Gi(t < 0) = 0, the Al - 2'th 
derivative of Gi(t) has a (5-function at t = 0, so that 



G,H ~ 1/w 



4Z-2 



i < N/2 



(23) 



for large u. 



DISCUSSION 



In this paper, we have given an exact linear response 
formula for the current in a wire in response to time- 
dependent temperatures applied at the boundaries. For 
a harmonic chain the full response function has been ana- 
lytically computed. We have presented numerical results 
for the frequency dependence of the current response in 
oscillator chains. For a diffusive system we find that the 



response differs from what is expected from a solution 
of the Fourier's equation with oscillating boundary tem- 
peratures. It is straightforward to generalize the deriva- 
tion to fluid systems, various stochastic and deterministic 
baths, and arbitrary system size L and spatial dimension 
d. This is discussed in detail in Ref. [13| for the u = 



As shown in Ref. 13j the zero frequency response can 
be expressed in terms of current auto-correlation func- 
tions, resulting in an expression similar to the standard 
Green-Kubo formula but without the thermodynamic 
limit being taken. If the integral of the auto-correlation 
function remains finite in the thermodynamic limit, the 
conductance is ~ 1/N for large N, and one can define 
an TV-independent conductivity in the same regime. The 
resultant expression matches the standard Green-Kubo 
formula, but with the thermodynamic limit taken after 
the range of the integral is taken to infinity. While it 
is plausible to assume that the order of limits commutes 
and 



lim — lim 

L— foo L to^oc 



Cjj(t)dt = lim lim — 



Cj,,(t)dt, 

this is by no means trivial: if different boundary con- 
ditions had been employed, with hard wall boundaries 
instead of heat baths, the left hand side of this equa- 
tion is zero but the right hand side is not. If the left 
hand side (with heat bath boundary conditions) diverges 
in the thermodynamic limit, as for integrable systems or 
low dimensional momentum conserving systems, the con- 
ductivity also diverges, and one can only talk about the 
conductance or an L-dependent conductivity. 

At non-zero frequencies, the integral converges even 
when it does not at w = 0, and changing the order 
of limits is more benign. Unfortunately, as we have 
seen in this paper, the expression obtained for the finite 
frequency conductance involves the correlation function 
(i/+i,/( T )^/p(0))j which we are unable to convert into an 
auto-correlation function when lj ^ 0. The connection to 
proposed expressions for the finite frequency conductiv- 
ity H, [Io[ is not clear. 
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